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ABSTRACT 

The transition density p t and pressure P t at the inner edge separating the 
liquid core from the solid crust of neutron stars are systematically studied using 
a modified Gogny (MDI) and 51 popular Skyrme interactions within well estab- 
lished dynamical and thermodynamical methods. First of all, it is shown that 
the widely used parabolic approximation to the full Equation of State (EOS) 
of isospin asymmetric nuclear matter may lead to huge errors in estimating the 
transition density and pressure, especially for stiffer symmetry energy functionals 
E sym (p), compared to calculations using the full EOS within both the dynamical 
and thermodynamical methods mainly because of the energy curvatures involved. 
Thus, fine details of the EOS of asymmetric nuclear matter are important for 
locating accurately the inner edge of the neutron star crust. Secondly, the tran- 
sition density and pressure decrease roughly linearly with the increasing slope 
parameter L of the E sym (p) at normal nuclear matter density using the full EOS 
within both the dynamical and thermodynamical methods. It is also shown that 
the thickness, fractional mass and moment of inertia of neutron star crust are 
all very sensitive to the parameter L through the transition density p t whether 
one uses the full EOS or its parabolic approximation. Moreover, it is shown that 
the E sym (p) constrained in the same sub-saturation density range as the neutron 
star crust by the isospin diffusion data in heavy-ion collisions at intermediate 
energies limits the transition density and pressure to 0.040 fm~ 3 < p t < 0.065 
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fm -3 and 0.01 MeV/fm 3 < P t < 0.26 MeV/fm 3 , respectively. These constrained 
values for the transition density and pressure are significantly lower than their 
fiducial values currently used in the literature. Furthermore, the mass-radius 
relation and several other properties closely related to the neutron star crust are 
studied by using the MDI interaction. It is found that the newly constrained pt 
and Pt together with the earlier estimate of AI/I > 0.014 for the crustal fraction 
of the moment of inertia of the Vela pulsar impose a more stringent constraint of 
R > 4.7 + 4.OM/M km for the radius R and mass M of neutron stars compared 
to previous studies in the literature. 

Subject headings: transition density — symmetry energy — stars: neutron - 
stars: crust 



1. Introduction 



Neutron stars are among the most mysterious objects in the Universe. They are nat- 
ural testing grounds of our knowledge about the Equation of State (EOS) of neutron- 



tron stars ( 


Lattimer & Prakash 


2007 


Ovamatsu & Iida 


2007; 


Douchin & Haense 




2000; 


Horowitz & Piekarewicz 


2001; 


Kubis 


2007; 


Ducoin et al. 


2007 


Rabhi et al. 


2008 


). Neu- 



tron stars are expected to have a solid inner crust which is belie yed to play an impo r- 



tant role in understanding a number of as 


jtrophysical observations ( 


Bavm e 


t al. 


1971a 


b; 




Pethic 


t & Ravenhall 


1995; 


Pethick et al. 


1995: 


Lattimer & Prakash 


2000: 


Steiner et al. 




2005a: 


Lattimer & Prakash 


2007 


Chamel & Haensel 1 


2008 


), such as, pulsar glitches 
m following x-ray bursts 


(Link et al. 


1999 


), quasi-period] 


c oscillations observed in x-ray emissic 


on neu- 



tron star (IDuncan 1 119981 ). the cooling observed over the first several years followin g: su 



per bursts from neu tron stars or giant flares from magnetars (Rutledge et al. II2OO6I ). and 



neutrino opacities (jHorowitz et al. 1 120041 ; iBurrows et al. 1 120061 ) . The solid inner crust of 



a neutron star comprises the region between the density p out where neutrons drip out of 
nuclei and the density p t where the transition to the homogeneous nucleonic matter occurs. 
e the p n ut is relatively well determ ined to be p out ~ 4 x 10 11 g/cm 3 (IRuster et al. 



Whi] 


c 


2006 





Hempel fc Schaffner-Bielich 1 120081 ) . the transition density p t is still very uncertain 



( ILattimer Sz Prakash Il2000l . 120071 ). This is largely due to our poor knowledge about the EOS 
of neutron-rich nuclear matter, especially the density dependence of t he nu clear symmetry 
energy E sym (p) at sub-saturation densities (ILattimer &: Prakash Il2000l . 120071 ). Consequently, 
our ability of understanding accurately many impo rtant properties of neutron stars has been 
hampered (ILattimer fc Prakash lliool . I2OO0I . l2007h . 
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The EOS of neutron-rich nuclear matter also plays an important role in heavy-ion colli- 
sions especially those induced by neutron-rich radioactive beams in terrestrial laboratories. 
While heavy-ion collisions are not expected to create the same matter and conditions as in 
neutron stars, the same elementary nuclear interactions are at work in the two cases. Thus, 
it is important to examine ramifications of conclusions regrading the EOS extracted from one 
field in the other one. Significant progress has been made recently in constr aining: the EOS 
of neutron-rich nuclear matter using heavy-ion experiments (See, e.g., ref. (ILi et al. 1120081 ) 
for the most recent review). In particular, compar ed to the existing model predictions in 
the literature the analyses of isospin diffusion data flTsang et al. 1120041 ; IChen et al. Il2005al ; 



Li fc Chen 1120051 ) in heavy-ion collisions have constrained relatively tightly the E sym (p) in 



exactly the same sub-saturation density region around the expected inner edge of neut ron 



star crust. Moreover, concl usions from analyz ing; some recent data (jShetty et al. 1 120071 ) of 



the isoscaling phenomenon (ITsang et al. 



on the thickness of neutron-skin in 208 Pb ( iSteiner fc Li Il2005bl ; ILi fc Chen 112005k I Chen et al. 



20011) in heavv-ion collisions and the available data 



2005bl ) are consistent with the E sym (p) constrained by the isospin diffusion data. Further- 
more, the lower bound of the experimentally constrained E sym (p) is consistent with the 
Relativistic Mean Field model prediction using the FSUGold interaction that can reproduce 
not only saturation properties of nuclear matter bu t also structure properties and giant 
resonances of many finite nuclei (IPiekraewicz 1 120071 ). While some model dependence and 
uncertainties still exist in the analyses of the above mentioned experiments and calculations, 
an overlapping area of the e xtracted E xym (p) from several analyses has appeared in the 



sub-saturation density region (ITsang et al. 1120081 ; iLynch et al. II2009I ). On the other hand 



extremely impressive progress has also been made in astrophysical observations relevant for 
constraining the EOS of nuclear matter. To our best knowledge, nevertheless, mainly be- 
cause of the low precision associated with the current measurements of neutron star radii, 
a non-controversial conclusion on the EOS and the density dependence of symmetry energy 
has yet to come. More accurate observations of neutron stars properties, especially their 
radii, with advanced x-ray satellites and other observatories, will hopefully enable us to con- 
strain stringently the EOS of neutron-rich matter in the near future. A direct cross-check 
on the EOS extracted independently from heavy-ion reactions and neutron star observations 
will then be possible. In the meantime, examinations of astrophysical implications of the 
EOS constrained by heavy-ion reactions are useful. At the WCI3 meeting in 2005, Horowitz 
suggested the heavy-ion physics community to investigate whether one can use the infor- 
matio n from heavy-io n collisions to constrain the core-crust transition density in neutron 
stars ([Horowitz 1 120051 1 . It is thus interesting to investigate timely how the behaviors of the 
E sym (p) constrained at sub-saturation densities by heavy-ion experiment s may help limit t he 
transition density p t and pressure P t at the inner edge of neutron stars (IXu et al. Il2008bl ). 
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To our best knowledge, all existing studies indicate consistently that the transition den- 



sity i s very sensitive to the density dependence of the nuclear symmetry e nergy (ILattimer fc Prakash 



20071 : lOvamatsu fc Iida 1 120071 : iDouchin fc Haensel 11200(1 Kubis 1120071 ). Very often, the so- 



called parabolic approximation (PA) to the EOS of isospin asymmetric nuclear matter is 
used. While the PA is mathematically valid only at small isospin asymmetries, interestingly, 



ing most models and interactions, see, e.g., refs. ( 


Bombaci & Lombarc 


o 


1991; 


Chen et al. 


2001; 


Zuo et al. 


2003; 


Xu et al. 


2007a: 


van Dalena et al. 


2007; 


Mou 


2007 


). Nevertheless, 



since the npe matter in the crust at /5-equilibrium is highly neutron rich and the determina- 
tion of the transition density depends on the second order derivatives of the energy density, 
the fine details of the EO S can influence the transition density significantly as first pointed 
out by Arponen in 1972 (lArpoen Ill972l ). It is thus interesting and necessary to compare 
calculations using both the full EOS and its parabolic approximation. Indeed, we found 
that the PA leads to significantly different transition density and pressure compared to the 
calculations using the full EOS. It should be mentioned that the PA may also significantly 
modify the proton fraction in /^-equilibrium neutron-star matter and the cr itical density for 



the direct Urea pro cess which can lead to faster cooling of neutron stars (IZhang fc Chen 



200ll ; ISteiner II2006I ). To investigate effects of nuclear interactions we use a modified Gogny 
(MDI) and 51 Skyrme interactions widely used in the literature. The same MDI interac- 
tion ha s been used in extracting the E svm (p ) from heavy-ion reactions within a transport 



model flChen et al. Il2005al; iLi fc Chen 1 120051 ) . Using the E sym (p) constrained by the isospin 
diffusion data fjTsang et al. 1120041 ). we can put a constraint on the transition density and 
pressure, respectively. We will then examine the implications of these constraints on the 
mass-radius correlation and the crustal fraction of the moment of inertia of neutron stars. 

This paper is organized as follows. In Section [2] we briefly review the dynamical and 
thermodynamical methods widely used for locating the inner edge of neutron star crust, and 
derive their relationship analytically. In Section [3] we summarize the EOS and symmetry 
energy obtained using the MDI interaction and 51 Skyrme interactions within the Hartree- 
Fock approach. We also examine the associated proton fraction and several thermodynamical 
properties including the energy density, pressure and the speed of sound in neutron star 
matter at /5-equilibrium. The general formalisms for describing the structure of neutron 
stars are outlined in Section HI We thus present the results of our calculations and discuss 
several important issues regarding the transition density and the structure of neutron stars 
in Section [51 A summary is given in Section [6j 



- 5 - 



2. Methods for locating the inner edge of neutron star crust 



The inner edge of neutron star crust corresponds to the phase transition from the 
homogeneous matter at high densities to the inhomogeneous matter at low densities. In 
principle, the transition density p t can be obtained by comparing a detailed model of 
the nonuniform solid crust to the uniform liquid core in the neutron star. While this is 
practically very difficult s ince the inner crust may have a very complex structure, usually 



known as "nuclear pasta" (Ravenhall et al. 111983 



2008a 



Ovamatsu Hl993l ; lHorowitz et al. 



200 



4; 



Steiner 



b|), it can be explored within several a 



...Hashimoto et al. 111984 
20081 : iGoeelein et al. ~ 



Lorenz et al. 



2008 



1993 



Avancini et al. 



aproaches including the molecular dynamics 



simul ations (IWatanabe 1 120051 ; iHorowitz 1 120061 ) and the 3D Hartree-Fock model ((Newton 



20091 ). Furthermore, the core-crust transition is thought to be a very weak first-order phase 
trans ition and model calculations lead to very smal l dens i ty discontinuities at the transi- 



tion ((Petliicket^ [1J9J; [Dou^^ Alterna- 



tively, a well established approach for estimating the p t is to search for the density at which 
the uniform liquid first becomes unstable against small-amplitude density fluctuations, indi- 
cating the start of forming nuclear clusters. Although some quantum effects such as the shell 
effects in more microscopic methods may influence the core-crust transition density, this 
approach has been shown to produce a very small error for the actual core-crust transition 
densi ty and it would yield the exact transition density for a second-order p hase transi- 



tion flPethick et al. Ill995t iDouchin &: Haensel 11200(1 l200ll ; ICarriere et al. 



there are several such methods, such as, the dynamical method (|Bavm et al 



Pethick fc Ravenha 



2007 



2007 



Ducoin et al. 



Worley et al. 



1995; Pethick et al. Il995: Douchin fc Haensel 



2007), the thermodynamical method (Kubis 



2008a 



Kubis 



2007 



Lattimer fc Prakash 



Approximation (RPA) (IHorowitz fc Piekarewicz 



2001 



2003h. 



2000 



2007 

T 



'resent 



1971a 



Ovamatsu fc Iida 



Lattimer fc Prakash 



2007) and the Random Phase 



Carriere et al. II2003I ). In the present 



work, we use both the dynamical and thermodynamical methods. 

In the following, we will first review briefly the dynamical method and the thermody- 
namical method, separately. While they are both well established and applied extensively 
in studying not only the core-crust transition in neutron stars but also the liquid-gas phase 
transition in asymmetric nuclear matter, somewhat different results are often obtained. It 
is thus necessary to study in detail the differences and relations between them. We shall 
first show analytically that the thermodynamical method corresponds to the long-wavelength 
limit of the dynamical method when the Coulomb interaction is neglected, and then compare 
numerically their predictions. 
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2.1. The dynamical method 



To describe small d ensity fluctuations in the npe matter, one can write the dens ity of 
particle q G {n,p,e} as fjBaym et al. Ill971at iPethick et al. Ill995l ; iDucoin et al. 1120071 ) 



Pq = Pg + 5pq- 

The density variation can be decoupled into plane-waves 



A q e ik - r + c.c, 



(2) 



of wave vector k and amplitude A q . This kind of density variation occurs when a momentum 
k is transferred to the particle system, e.g., through collisions and the "dynamical method" 
is named after this. It has been shown that the variation of the free energy dens ity generated 
by the density fluctuation of amplitude A = (A n , A p , A e ) can be written as (IBaym et al. 



1971al : IPethick et al. Ill995t IDucoin et al. 1120071 ) 



where 



5f = A*C f A, 



dfi n /dp n dfi n /dp p 
dn P /dp n dp p /dp p 

Ofle/dpe 



(3) 



D n n D n p 
Dp n Dpp 



Aire' 












(4) 



is the free-energy curvature matrix. The instability region of the npe matter can be located by 
examining when the convexity of the free-energy curvature matrix is violated. The convexity 
of the matrix requires that 



C{ x > or C( 2 > 0, 



nf 



rJ 



>o, 



rif r<f nf 

°11 °12 °13 

rif rif rif 

21 °22 °23 

rif rif rif 

u 31 ^32 ^33 



> 0. 



(5) 



Here C( 3 is always positive so we do not take it into consideration. If the system stays 
stable, the convexity of the matrix C* should be retained for all values of k. The first 
term in the right hand of Eq. (j3J) is the bulk term, which just defines the stability condition 
of the nuclear matter part as will be shown later. The second term in the right hand of 
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descri bes the contribution o f the density gradient. For the Skyrme-Hartree-Fock 



Eq. I 

(SHF) model fjChabanat et al. 



D, 



D 



D 



Tip 



D,, 



19971 ) one has 
3 

16 
1 

16 



pp 



pn 



[tl(l-Xi)-t 2 (l + X 2 )], 

[3t 1 (2 + x 1 )-t 2 (2 + x 2 )] 



(6) 



(7) 



in terms of the standard Skyrme interaction parameters xi,x 2 ,ii and t 2 - The MDI inter- 
action, however, does not have a gradient term. To remedy this draw back we set D pp = 



D 



up 



132 MeV-fm 5 as used in the work by Oyamatsu et al. (jOyamatsu fc Iida 



20071 ) when we apply the MDI interaction. This choice is quite consistent with the empirical 



values from the SHF calculations. We note here that the averaged value of D pp = D nn and 
D np = D pn is, respectively, 140.9 and 118.8 MeV-fm 5 for the 51 Skyrme forces we will use 
in the following. Furthermore, as we will show later, the transition density and pressure 
are rather insensitive to the variation of D pp = D nn and D np = D pn . The last term in the 
right hand of Eq. (j3J) is the Coulomb term, which is generated by the Coulomb interactions 
of electrons and protons. It should be noted that additional /c-dependent terms due to the 
finite range of the MDI interaction via exchange terms as well as the Coulomb exchange 
terms are neglected in Eq. (j4j). As we will show later, the bulk term dominates the result 
and the density gradient term and Coulomb term are not important for the determination of 
the transition density and the associated transition pressure. The density gradient term and 
Coulomb term usually make the system slightly more stable and thus reduce correspondingly 
the region of instability. 

For small density fluctuations, to guarantee the con yexity of the curvature matrix it is 
suffic ient for the last determinant in Eq. ([5]) to be positive (IBaym et al. Ill971bl ; iPethick et al. 
1995h . i.e., 

V + (3k 2 + >0, 



where 



V 





Vdyn(k) PS 

dp p {dp n /dp v 



k 2 + B 



TF 



dpp 
D 



dji n /dp n 



pp 



2D np ( + D nn ( 2 , ( 



dp p /dp n 

dp n /dp n 



k 



2 

TF 



Aire 2 

dpe/pe 



(9) 
(10) 

(11) 



dp p \dp 



In the above expressions, we used the relation §^ = |^ following 

1 ' ap p ap n ° dpp 

ef - (^") = ^ W ^ n e t» ern g the energy density of the npe matter. Meanwhile, is 
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assumed to be positive. If we have -J^- < but -Jr 2 - > we can change the form of the 
equations correspondingly. In this form, it is clear that the density gradient and the Coulomb 
term clearly contribute positively to the Vd yn (k). They thus help to m ake the system more 



stable. At k=\(^f -) 1 ' 2 - k 2 TF } 1 / 2 , the V r 



Pethick et al. 1119951 ) 



dyn 



(k) has a minimal value of flBaym et al. Ill971b 



V dyn = V + 2(4vre 2 /3) 1 / 2 - (3k 2 TF . (12) 
Then the density at which Eq. (|12p becomes zero determines the instability boundary. 



2.2. The thermodynamical method 



The thermodynamical metho d requires the system to obey the stability condition (IKubis 
20071 : lLattimer fc Prakash~ll2007h 

(13) 



dv , 



dp, 
dq c 



> 0, 



or the system will be unstable against small density fluctuations. These conditi ons are 



(14) 



equiv- 



alent to requiring the cony exity of the energy per particle in the single phase (IKubis 1 120071 : 
Lattimer fc Prakash 1 120071 ) by ignoring the finite size effects due to surface and Coulomb 
energies as shown in the following. Here the P = P b + P e is the total pressure of the npe 
system with the contributions Pb and P e from baryons and electrons, respectively. The v and 
q c are the volume and charge per baryon number. The \x is the chemical potential defined as 



/* (J"n ftp • 



(15) 



In fact, Eq. (flBj) is simply the well-known mechanical stability condition of the system at 
a fixed jj,. It ensures that any local density fluctuation will not diverge. On the other 
hand, Eq. (IT41) is the charge or chemical stability condition of the system at a fixed density. 
It means that any local charge variation violating the charge neutrality condition will not 
diverge. If the /3-equilibrium condition is satisfied, namely \x = ft e , the electron contribution 
to the pressure P e is only a function of the chemical potential /i, and in this case one can 
rewrite Eq. (TTBl as 

'op; 



dv 



> 0. 



(16) 



By using the relation dEb ^' Xp ^ = — ^ one can get f Kubis 2007 ) 



dv 



2p 3 



dE b (p,x p 
dp 



fi 2 E b (p, x f 
dp 2 



d 2 E b (p,x r 
dpdxp 



2 , d 2 E b ( y p,x p ) 
dx 2 



(17) 



- 9 - 



d 2 E b (p,x p ) dp 



dx 2 



+ 



where g c 



p e l p. The p = 1/v is the baryon density and the E b (p,x p ) is the energy per 



baryon. Within the free Fermi gas model, the density of electrons p e is uniquely determined 
by the electron chemical potential Up- Then the thermodynamic al relations Eq. (TTBl and 
Eq. (P3J) are identical to jLattimer fc Prakash 1120071 : kubis 1 120071 ) 



m 

dv 



2 P 



dE b (p,x f 
dp 



2 d 2 E b (p,x p ) 
dp 2 



d 2 E b (p,x p ) \ d 2 E b (p,x p 



dpdxp 



P I- 



dx 2 



= / d 2 E b (p,x p ) | > 

\ dp J v dx 2 ir 2 h 3 p 



> 0,(19) 



(20) 



respectively. The second inequality is usually valid. Thus, the following condition from the 
first one 



T/ o dE b (p,x p ) 2 d 2 E b (p,x p ) 

Vther = 7{ h P 



d 2 E b (p,x PJ ^, 



dp dp 2 \ dpdxp 

determines the thermodynamical instability region. 



d 2 E b (p,x p ) 
dx 2 



(21) 



Within the parabolic approximation neglecting higher order terms of isospin asymmetry 
5 = 1 — 2x p , the EOS of asymmetric nuclear matter is 



E b (p, S) = E (p) + E sym {p)8 2 



(22) 



where Eq(p) is the energy per nucleon of s ymmetric nuclear matter. Then Eq. ([2"T]) can be 
reexpressed as ( jLattimer &: Prakash 1 120071 ) 



„ d 2 E dE 

ther — P ~T~T + l P~ 



yPA _ „2 



dp 2 



dp 



2X r 



p 



d 2 E dE 



dp 2 



dp 



sym 



dE 



sym 



dp 



(23) 



2.3. The relationship between the dynamical and thermodynamical methods 

The Eq. ( Tl2l) and Eq. ( 12T1) together with the relationship between the density p and the 
proton fraction x p required by the /3-equilibrium and the charge neutrality conditions will 
then determine respectively the dynamical and the thermodynamical core-crust transition 
density in neutron stars. These two methods together with various EOS's have been widely 
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used in the literature while their relationship is still unclear. Therefore, it would be inter- 
esting to first obtain some analytical insights on their relationship before comparing their 
numerical predictions. 

In the following, we first analyze the instability of asymmetric nuclear matter without 
considering the /3-equilibrium and the charge neutrality conditions. As used in the previous 
subsection and the literature, the stability condition is often expressed using the p and x p 
within the dynamical method while the p n — (1 — x p )p and p p = x p p within the thermody- 
namical one, respectively. Thus the following simple thermodynamical relations are useful 
for understanding the relationship between the two methods 



dEb 

dx„ 



- fi flp P"ni 



m 

dp 

m 

dx„ 



x p )p 
x p )p 



dp 

dfi n 
dx n 



x p p 



Xpp 



d/ip 

dp 

dfip 

dx n 



(24) 

(25) 
(26) 



where the pressure of baryons is Pb = fi n p n + PpPp 
energy of baryons can be expressed as 



dE b 


P h 




dp 






d 2 E b 


= d_(Pj> 




dp 2 


dp \p 2 


) 




2P b 


1 




P 3 






2P b 


i 




P 3 


7 



E b p. In this way, the derivatives of the 



(27) 



(1 - x p)P-^J + x pP~^7 



dp 

(1 - Xpfp^^ + x p (l - Xp)p 



dpr, 



+ ~~n 



X , 



dpr, 

x v)p^t + X P P- 



dpn 

dp p 



d 2 E b 
dx 2 



d 2 E b 
dpdxp 



dp, dp, p 



dp n 

dx 



P 



dfip _ dfip _ djin 
dpp dp n dpp 

dfi dfip dfi n 
dp 



dpp 



dfi n 
dpn 



{21 



(29) 



(1 



Xr 



dp 
dfip 
dpn 



X 



dp 
dp p 
dpp 



dfin dfi n 



dpn 



X li 



dpp 



(30) 
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As shown earlier, for nuclear matter without considering the Coulomb interaction one has 

dp p dp n ' 

From Eq. (12 7p . (I28p . (I29p . (I30p and (1311) . we can then obtain the following important equality 
2dE b d 2 E b + ^E 1 ^E 1 _ f d 2 E b \ 2 _dp n dp p fdpA 2 ^ 



p dp dx 2 dp 2 dx 2 \dpdx p J dp n dp p \dp p 

Therefore, for positive values of ^J?, the condition Eq. (fT9l) is simply equivalent to requir- 
ing a positive bulk term Vq in the Eq. (J4]). Since the transition density is usually in the 
sub-saturation density region where the -g^r > is valid for almost all model EOS's, the 
thermodynamical stability condition is thus simply the limit of the dynamical one as k — > 
(long-wavelength limit) when the Coulomb interaction is neglected. 



3. The EOS and symmetry energy with selected 51 Skyrme forces and a 

modified Gogny interaction 

In this section, we summarize the EOS and the corresponding s ymmetry energy o btained 
using the modified finite-range Gogny effective interaction (MDI) (IDas et al. 1 120031 1 and 51 
popular Skyrme forces within the Hartree-Fock approach. These results will be used later in 
our numerical calculations of the core-crust transition density and pressure. The MDI inter- 
action has been extensively used in our previous studies of heavy-ion collisions, the liquid-gas 



phase transition in asymmetric nuclear matter and several issues in astrophysics fjLi et al 



20081 ). The EOS's using various Skyrme forces are well known for thei r simple forms and suc- 
cessful descriptions of many interesting phenomenon, see, e.g., ref s. (IChabanat et al. 1119971 ; 
Steiner et al. Il2005al ; IStone et al. 1120031 ; IStone fc Reinhard 1 120071 ). A very useful feature of 
both the MDI and the Skyrme interaction is that analytical expressions for many interesting 
physical quantities in asymmetric nuclear matter at zero temperature can be obtained. 



3.1. The EOS and symmetry energy with selected 51 Skyrme interactions 

Within t he SHF approach the e nergy per nucleon for symmetric nuclear matter can be 
expressed as (IChabanat et al. II1997I ) 



Eo(p) 



3h 2 /3tt 



10m 



2\ 2/3 



J2/3 



i 3 „ /3tt' 

-t p H S — 

I H 80 V 2 



2\ 2/3 



.,5/3 



16 



hp' 



cr+l 



(33) 
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Table 1: Saturation density po (fm -3 ), binding energy of symmetric nuclear matter Eq(pq) (MeV), 
incompressibility Kq (MeV), symmetry energy E sym (po) (MeV) as well as slope and curvature 
parameters of symmetry energy L (MeV) and K sym (MeV) at the saturation density. 



SHF 




Po 


£o(po) 


K 


E sy m{po) 


L 


K 


BSk3 


0. 


157 


-15.8 


234.8 


27.9 


6.8 


-306.9 


BSkl 


0, 


157 


-15.8 


231.3 


27.8 


7.2 


-281.8 


BSk2 


0, 


157 


-15.8 


233.7 


28.0 


8.0 


-297.0 


MSk7 


0, 


157 


-15.8 


231.2 


27.9 


9.4 


-274.6 


BSk4 


0, 


157 


-15.8 


236.8 


28.0 


12.5 


-265.9 


BSk8 


0, 


159 


-15.8 


230.3 


28.0 


14.9 


-220.9 


BSk6 


0, 


157 


-15.8 


229.1 


28.0 


16.8 


-215.2 


BSk7 


0, 


157 


-15.8 


229.3 


28.0 


18.0 


-209.4 


SKP 


0, 


163 


-16.0 


201.0 


30.0 


19.6 


-266.8 


BSk5 


0, 


157 


-15.8 


237.2 


28.7 


21.4 


-240.3 


SKXm 


0, 


159 


-16.0 


238.1 


31.2 


32.1 


-242.8 


RATP 


0, 


160 


-16.0 


239.4 


29.2 


32.4 


-191.2 


SKX 


0, 


155 


-16.1 


271.1 


31.1 


33.2 


-252.1 


SKXce 


0, 


155 


-15.9 


268.2 


30.1 


33.5 


-238.4 


BSkl5 


0, 


159 


-16.0 


241.6 


30.0 


33.6 


-194.3 


BSkl6 


0. 


159 


-16.1 


241.7 


30.0 


34.9 


-187.4 


BSklO 


0, 


159 


-15.9 


238.8 


30.0 


37.2 


-194.9 


SGII 


0, 


158 


-15.6 


214.7 


26.8 


37.6 


-145.9 


BSkl2 


0, 


159 


-15.9 


238.1 


30.0 


38.0 


-191.4 


BSkll 


0, 


159 


-15.9 


238.1 


30.0 


38.4 


-189.8 


BSkl3 


0, 


159 


-15.9 


238.1 


30.0 


38.8 


-187.9 


BSk9 


0, 


159 


-15.9 


231.4 


30.0 


39.9 


-145.3 


SLylO 


0. 


158 


-16.5 


237.8 


33.2 


40.8 


-148.0 


BSkl4 


0, 


159 


-15.9 


239.3 


30.0 


43.9 


-152.0 


SLy230a 


0, 


160 


-16.0 


229.9 


32.0 


44.3 


-98.2 


SKM* 


0, 


160 


-15.8 


216.6 


30.0 


45.8 


-155.9 



13 



Table 2: Continued with Table [T1 
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33 1 
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1 9fi Q 

-izo.y 


D1Y1V1 
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u 
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-lO.O 
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ZIO.O 


3D 7 
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ZLQ 3 
4y.o 
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00.4 


4Q 7 
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-1U.O 
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OO.O 


a 


1 9n 9 

-1ZU.Z 


oijyo 


n 
u 


1 fi3 
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-lO.O 


938 n 

ZOo.U 


33 fi 
OO.O 


^1 Q 
oi.y 


1 1 fi 3 
-1 lO.O 


SLy9 


0. 


153 


-16.4 


237.7 


33.2 


57.2 


-84.9 


SkI6 


0. 


159 


-15.9 


248.2 


29.9 


59.2 


-46.8 


SkI4 


0. 


160 


-15.9 


248.0 


29.5 


60.4 


-40.6 


SGI 


0. 


154 


-15.9 


261.8 


28.3 


63.9 


-52.0 


SKO* 


0. 


160 


-15.7 


222.1 


32.1 


69.7 


-77.5 


SkMP 


0. 


159 


-16.1 


238.5 


30.1 


70.7 


-51.4 


SKa 


0. 


155 


-16.0 


263.2 


32.9 


74.6 


-78.5 


SKO 


0. 


160 


-15.8 


222.8 


32.0 


79.5 


-42.3 


R a 


0. 


158 


-15.6 


237.4 


30.6 


85.7 


-9.1 


SKT4 


0. 


157 


-15.5 


229.3 


34.8 


92.4 


-24.2 


G<j 


0. 


158 


-15.6 


237.2 


31.4 


94.0 


14.0 


SkI3 


0. 


158 


-16.0 


258.2 


34.8 


100.5 


73.0 


SkI2 


0. 


158 


-15.8 


240.9 


33.4 


104.3 


70.7 


SkI5 


0. 


156 


-15.8 


255.8 


36.6 


129.3 


159.6 
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wi th Q, s = 3ti + (5 + 4 x2)^2- For asymmetric nuclear matter, the energy per nucleon 
is ( IChabanat et al. II1997I ) 

ofc2 /o 2\2/3 1 

E b ( Pl 5orx p ) = ^(^M p 2/3 F 5/3 + -M2(x + 2)-(2xo + l)F 2 ] 

I o /o 2\2/3 

+ ^ +1 [2(x3 + 2)-(2x 3 + l)F 2 ] + A^j p5 /3 

x j^z! + 2) + t 2 (x 2 + 2)]F 5/3 + ^[t 2 (2x 2 + 1) - t!(2x! + 1)]F 8/3 | , 



(34) 



with 



F m (s) = \[(i+s) m +o--n 

F m (x P )=2 m - 1 [x p n + (l-x p ) m }. 

Within the parabolic approximation widely used in the literature, the symmetry energy is 
calculated from 

E sym (p) « E b (p, 5=1)- E b (p, 5 = 0). (35) 

But strictly speaking, the symmetry energy should be the coefficient of 5 2 in the Taylor 
expansion of E b (p, 5) in terms of 5, i.e., 

Enm(p) = i(^) (36) 



We notice here that the above two definitions for the symmetry energy would be the same 
should there be no higher order terms in 5 in the EOS of asymmetric nuclear matter (But 
it should be noted that the kinetic part of the EOS of asymmetric nuclear matter always 
contains higher order terms in 5). 

Thus, by definition of Eq. (13*51) . for Skyrme interactions, one has 
p / \ 1 (d 2 E, 

sym\P ) 



2 \ d5 2 



8=0 

2\ 2/3 



h 2 f^ 2 y /A 2/3 1 , 

(jy) Qs Vm p 5/3 -^t 3 (2x 3 + l)p° + \ (37) 



24 

where O sym = 3tiXx — 1 2 (4 + 5x 2 ). a, t ~ t 3 and x ~ x 3 are the Skyrme parameters. 
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As it has been used extensively by many authors, near the saturation density po the 
symmetry energy can be expanded as 



E S ym{p) ~ E sym (p ) 



L 

3" 



P- Po 
Po 



K. 



sym / P Po 



PO 



(3J 



where L and K sym are, respectively, the slope parameter and curvature parameter of the 
symmetry energy at p , i.e., 



L 



3p, 



K 



9Po 



dE sym (p) 
dp 

d 2 E sym (p) 



dp 2 



P=PQ1 



I p=po ■ 



(39) 
(40) 



The L and K sym can be used conveniently to characterize the density dependence of the 
symmetry energy around the saturation density p . In the present work we use 51 standard 
Skyrme forces with their saturation density and the symmetry energy satisfying 0.140 fm~ 3 < 
Po < 0.165 fm -3 and 26 MeV< E sym (po) < 37 MeV, respectively. Some Skyrme forces with 
very small or negative L values are not considered here as they generally predict bound pure 
neutron matter at sub-saturation densities and are not suitable for the descriptio n of neutron- 



rich e nvironments like neutron star crusts as discussed in detail by Stone et al. (IStone et al 



20031 ). In addition, we have not included the Skyrme forces predicting values for the incom- 



values of the parameters for these 51 Skyrme forces can be found in refs. (I 


3rack et al. 


1985; 


Friedrich & Reinhard 1986: 


Brown 1998; Chabanat et al. 


1997: Stone et al. 


2003 


Stone & Reinhard 1 


2007: Chen et al. 


2005bl: Steiner & Li 1 2005b: 


Samvn et al. 


2002. 


20031: 


Gorie 


lv et al. 12003 


Samvn et al. 


200 


4;lGorielv et al. 


2005: Samvn et al. 2005 


Gorielv et al. 


2006. 


2007: Chamel et al. 2008a: 


Gorielv & Pearson 


20081). The selected ranges of pn and 



E sym (p ) are consistent with their empirical values inferred from nuclear laboratory data. 
The detailed properties of these forces at p are summarized in Tables [1] and [2] in the order 
of rising values of L. 



3.2. The EOS and symmetry energy with the modified Gogny interaction MDI 



For the MDI interaction based on the Hartree-Fock calculat ion using the Go gny inter- 
action, the baryon potential energy density can be expressed as (IDas et al. 1120031 ) 



V(p,6) 



nPp 



Po 



+ ^-^(Pn + P P ) + 



2p 



*+l Po 



+ — c tt i 
po f-f 




d 3 pd 3 p 



1 + (p-^) 2 /A 2 



(41) 
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We notice here that the above is a natural extension to isospin asymmetric ca se of the cor- 



responding; potential energy density for symmetri c nuclear matter giv en in refs. flGale et al 



19871 : IPrakash et al. I Il988l : IWelke et al. I Il988l : iGale et al. I Il990l ). 



BGBD (Bombaci-Gale-Bertsch-Das Gupta) potential energy density (IBombaci 



It is similar to the 
The 



gas phase transitions in neutron-ri ch matter (Xu et al. 112007a 



2001 



MDI interaction has been used extensively in stud ying heavy-ion reactions (ILi et a 
liquid 



Li et al. 



Li & Steiner 


2006 


Krastev & Li 


20C 


7 


Krastev et al. 


(Krastev et al. 


2008b: 


Worlev et al. 


2008b 


) of neu- 



2007 



2008). 



Xu et al. 



tron stars. 

In the mean field approximation, Eq. fill) leads to the following single particle potential 
for a nucleon with momentum p and isospin r, i.e., 



U(p,S,p,r) = A^x^ + Mx^ + B^Yil-xd 2 ) 

Po Po Po 



Pr 



STX- 



(7 



+ 1 (K ' 



+ 



2a,. f d3p) Uim + 20^ r d3p , f. T (f, f) 



Po 



l + (p-f) 2 /A 2 



Po 



1 + (p-f) 2 /A- 



, (42) 



In the abo ye the isospin r = 1 /2 (—1/2) for neutrons (protons). The coefficients A u (x) and 
A{x) are jGhen et al. Il2005al ) 



AJx) 



-95.98 -x- 



2B 

'<7+T 



Mx) 



-120.57 + x 



2B 
a+1 



(43) 



The values of the parameters are a = 4/3, B = 106.35 MeV, C T)T = —11.70 Me V, C 



103.40 MeV and A = p° which is the Fermi momentum of nuclear matter at po ( )Das et al. 



20031 ). The parameter x was introduced to mimic various E sym (p) predicted by different mi- 
croscopic many-body theories. By adjusting the x parameter, the E sym (p) is varied without 
changing any property of symmetric nuclear matter and the symmetry energy at saturation 
density as the x-dependent A u (x) and A[(x) are automatically adjusted accordingly. We 
note especially that the symmetry energy at normal density E sym (p Q ) is fixed independent 
of the x parameter. Using the definition in Eq. ( l36l) . E sym (p ) = 30.54 MeV at po = 0.16 
fm -3 while its value is 31.6 MeV within the parabolic approximation of Eq. (1331) . 

At zero temperature the phase space distribution function can be writte n as .f T (r, p) 



^■6(p/(r) — p), and all the integrals expressions can be calculated analytically (IWelke et al. 
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Fig. 1. — (Color online) The density dependence of the nuclear sym metry energy for diffe rent 
values of the parameter x in the MDI interaction. Taken from ref. (IChen et al. Il2005al ) 



1988; 


Das et al. 


2003; 


Chen et al. 


2007) 



p 2 f (r) + A 2 -^ in [p + J5/ ( r )] 2 + A 2 



2pA 



+ ^-2{arctan P + ^ (r) 



[p-p f (r)} 2 + A 2 

p-Pf(r) 



A 



A 



arctan : 



A 



}]• 



(44) 
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I (f) A2 ^/( r ^( r/ )[ 3 ^^) + - A 2 



+ 4A[(p^(r) - ?>/(t )) arctan 



A 



P/( r ) + P/( r ')- 



A 



x In 



(P/( r ) + P% T )) arctan 

^[A 4 + 6A 2 ( P J (r) + p 2 (r')) - 3(pJ(r) - pj(r')) 21 
(P/(r)+P/(r')) 2 + A 2 



(P/M-P/(r')) 2 + A 2 



}■ 



(45) 



The kinetic energy is 



E k (p,6) 



d 3 p 



An 



5mh*p {Pfn + Pfp) 



(46) 



where Pf n {p) = ^(37r 2 p n ( p )) 1 / 3 is the Fermi momentum of neutrons (protons). Then, the total 
energy per baryon for cold asymmetric nuclear matter is 



E b (p,5) = ^ P ^- + E k (p,S). 
P 



(47) 



By setting p n = p p = | and pf n = pf p = pf we thus obtain the following EOS of cold 
symmetric nuclear matter 



B 



P 



a+1 \po 



3poP 



(Q + C u ) 



x 



2p/ 



P/(6pJ - A ) - 8ApJ arctan + -(A + 12A 2 p^) In 



Ap 2 f + A 5 
A^ 



(48) 



We stress here that since the Ai(x) + A u (x) is a constant of —216.55 MeV according to 
Eq. (j43|) . the Eo(p) is independent of the parameter x as expected. The symmetry energy 
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by definition is 

E S ym (p) 



1 (d 2 E 



2 V d5 2 



8=0 



57r .5 , P i A A Bx ( p 



+ 



+ 



9mh 3 p 1 4pi 



a+l \po 



9p p \ h 3 J 
C u 



Ap 2 f + A 2 



9poP 



(49) 



where p/ = h^Ti 2 ^) 1 ^ 3 is the Fermi momentum for symmetric nuclear matter. We note here 



that since the A\(x) — A u {x) = —24.59 + ABx/(a + 1) according to Eq. 



the E sym (p) 



depends linearly on the parameter x at a given density except po where the symmetry energy 
is fixed by construction. As shown in Fig.[U adjusting the parameter x in the MDI interaction 
leads to a broad range of the density dependence of the nuclear symmetry energy, similar to 
those predicted by various microscopic and/or phenomenological many-body theories. 



3.3. Thermodynamics quantities in neutron stars at /3-equilibrium with the 

MDI interaction 

Since we are going to examine astrophysical implications of the symmetry energy con- 
strained by heavy-ion reactions obtained from transport model analyses using the MDI 
interaction, it is useful to first study several key thermodynamical quantities in neutron star 
matter at /5-equilibrium with charge neutrality. It is also necessary to examine the causality 
with the MDI interaction. 

It is well known that for the npefi matter the /3-equilibrium condition is 

Pn - Pp = Pe = AV (50) 

The appearance of muons requires a sufficiently high chemical potential of electrons, i.e. 
ix e > m^, where is the mass of muons. Eq. (150"]) together with the charge neutral 
condition 

pp = pe+Pv (51) 

determines the proton fraction x p as a function of baryon density in the neutron star matter. 

To calculate the core-crust transition density p t , we only need to deal with the npe 
matter since muons will normally not appear as the electron chemical potential p e is not 
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high enough near p t unless one uses an extremely soft symmetry energy. For the npe matter 
at /^-equilibrium, one has 

fin - A*p = A*e- (52) 

Then, this identity together with the charge neutral condition 

P P = Pe (53) 

gives the corresponding function of baryon density. 

If analytical expressions of the EOS for asymmetric nuclear matter are known completely 
as given earlier for the Skyrme and MDI interactions, the exact /3-equilibrium condition of 
Eq. (|21|) can be used. However, often this is impossible with many interactions within various 
models. Instead, the parabolic approximation of the EOS is usually used. In this case one 
has 

A* e = fJ>n ~ 4(1 - 2x p )E sym = A5E sym . (54) 

Using both the full EOS and its parabolic approximation of the MDI interaction with 
x = and x — — 1, we have calculated the proton fraction function of density from 

to 1.6 fm -3 . The specific values of the x pa rameter chosen here are consistent with the 



constraints extracted from heavy-ion reactions fjLi et al. II2008I ). The calculated values of x p 
are shown in the panel (c) of Fig. [2J Compared to the results with x — — 1, the x p with 
x = is larger below the saturation density and smaller at higher densities. The difference 
between calculations using the full EOS and its parabolic approximation is only visible at 
low densities for the soft symmetry energy with x = 0. 

We now examine several thermodynamical quantities for the npep matter at /3-equilibrium. 
The total energy density e(p, 8) consists of three parts: the baryon energy density e&(p, <5), 
the electron energy density e e (p, 5) and the muon energy density e M (p, 5) 

e(p, 5) = e b (p, 5) + e e (p, 5) + e M (p, 5), (55) 

where 

e b (p,8) = pE b (p,5) + pm (56) 

with m being the baryon mass and p the total baryon density. The energy density of leptons 
€i( p,5) is calculated using the no n-interacting Fermi gas model and it can be expressed 



as (jOppenheimer fc Volkoff Ill939l ) 



e l (p 7 5) = ri<f>(t), (57) 

with 

mic 2 
V = 8^' 
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Fig. 2. — (Color online) The density dependence of the energy density (a), the pressure (b), the 
proton fraction (c) and the sound velocity (d) for MDI interaction with x = and x = — 1 for the 
npefi matter at /^-equilibrium. The results from the full EOS and its parabolic approximation (PA) 
are compared. 



and 



A 



_h_ 

mic 



t = A(3vr>) 



<f>(t) = tVT+P(l + 2t 2 ) -ln(t + v / TTt2), 

where mi and pi are the mass and number density of leptons. 

Correspondingly, the total pressure P(p, 5) consists of the contributions from baryons, 
electrons and muons, i.e., 



P(p, 5) = P b (p, 5) + P e (p, 5) + P,(p, 5), 



where 



P b (p, 5) = p! n p n + p' p p p - e b (p, 5), 
and here the chemical potentials should include the rest mass 

p' n = Pn + m, p' p = p p + m. 



(58) 
(59) 

(60) 
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The pressure of leptons is written as 

p i(p>5) = fiipi - ei(p,5), (61) 

where the chemical potential is 



M = yPfi + m h ( 62 ) 
which is fully determined by the lepton density from 

Pfl = H(37r 2 Pl f 3 . (63) 

Then, in this framework the thermodynamical consistency 

P = ^ (64) 

is satisfied. 

The exact expressions given above can be carried out using the full EOS. While in some 
cases, the parabolic approximation is used. Instead of Eq. (j55l) and Eq. (j59]h within the 
parabolic approximation one has 

e b (p, 5) = p[E (p) + E sym (p)5 2 } + pm, (65) 

P b (p,5) = p 2 (E' (p) + E' sym (p)5 2 ). (66) 

This approximation still satisfy the thermodynamical consistency (Eq. (164)) ). 

The density dependence of the total energy density and pressure are shown in Panel (a) 
and Panel (b) of Fig. [21 respectively. The difference between calculations using the full EOS 
and its PA is essentially invisible. The stiffer (e.g., x=-l) the symmetry energy is, the larger 
the total energy and pressure are as one expects. 

The causality requires that the sound speed s in nuclear matter remains smaller than 
the speed of light in vacuum c, i.e., 



In Panel (d) of Fig. [2] we examine the speed of sound for the MDI interaction with x = 
and x = — 1. It is seen that the causality is satisfied in the whole density range considered. 
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4. Key equations for describing the structure of neutron stars 



F or completeness, we quote here from the general literature, see, e.g., ref. ( [Morrison et al. 
20041 ). some key equations to be used later in our studies of neutron star structure. For 
slowly-rotating neutron stars where the spherical symmetry is conserved approximately, the 
moment of inertia is 

' = <§>- = ^ < 68 > 

where Q is the angular velocity measured in a far-away inertial system and J is the angular 
momentum. In the slow-rotation limit in spherical polar coordinates the metric can be 
written as (G — c — 1) 

(ds) 2 = -e%dtf + {l-^)-\drf 

r 

- 2ur 2 sin 2 Odtdcj) + r 2 {d9 2 + sin 2 Odcj) 2 ), (69) 

where u{r) = ^ is the angular velocity of the local slow-rotation system (measured in a 
far-away inertial system), and m g (r) is the neutron star gravitational mass inside a radius r 

^M = A7rr 2 e(r), (70) 
dr 

with e(r) being the energy density. Defining 

ZJ = Q — u>, (71) 
then from the metric of the geometry outside a slow-rotation star one can get 

r -*jL[ r *j( r )^] + Ar- l ^-W{r) = 0, (72) 
dr L v ' dr i dr v ' v ' 

where 

j(r)=exp[--(A(r) + z/(r))], (73) 

with 

e A(r) = (1 _ ^hyl (74) 
r 

and 

dv 2 dP 

= (75) 

dr e + P dr' y ' 

The pressure P(r) is obtained from the famous Tolman-Oppenheimer-Volkoff (TOV) equa- 
tion, i.e., 

dP = m f + 4^P 

dr r(r — 2m g ) 



-24- 



For convenience, an additional function n(r) can be denned as 

,M = % (77) 
and the boundary conditions at the central of the star are 

uJ(0) = const, (78) 



and the constant is chosen so that 
Outside the star we should have 

and 

where 



r/(0) = 0. (79) 
2M 



r 

2J 



0) 



J = \R A n{R) (82) 

and the M and R are total gravitational mass and the total radius of the neutron star, 
respectively. To make the variables continuous at the surface of the star, we have the 
boundary conditions 

cJ{R) = tt- v(R)~, v{R) = ln(l - 2M/R). (83) 
3 

In this framework, after solving the differential equations Eq. flTUj) . (1721 . (176]) and (175]) . the 
mass, radius and moment of inertia can be calculated. Following the standard procedure, we 
integrate out the TOV and other differential equations equation starting from the center to 
the surface where the pressure vanishes, i.e., P(R) = 0. The latter defines the total radius R 
of the neutron star. Then the total gravitational mass of the neutron star is obtained from 
integrating Eq. (170]) as 

rR 

M = m g (R) = 4vr / drr 2 e(r). (84) 
Jo 

The total moment of inertia of the neutron star is obtained similarly. By integrating only 
to the transition density p t , one can obtain the radius and mass of the core. The thickness, 
mass and moment of inertia of the crust can be obtained from taking the differences between 
values for the whole and the core of neutron stars. 
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5. Results and discussions 



In the following, we present and discuss results of our calculations on the transition 
density and pressure at the inner edge and several global properties of neutron stars. Ap- 
plying formalisms outlined in the previous sections, we illustrate numerically and discuss 
several issues including (a) relationships among the mechanical, chemical and total instabil- 
ity boundaries in asymmetric nuclear matter, and their relevance for locating the core-crust 
transition density in the npe matter at /3-equilibrium; (b) the difference between the core- 
crust transition densities obtained using the dynamical and thermodynamical methods using 
the same interactions; (c) understanding the difference between the core-crust transition den- 
sities obtained with the full EOS and its parabolic approximation using the same methods 
and interactions; (d) the systematics of the transition density by varying the stiffness of 
the symmetry energy; (e) limits on the transition density using the symmetry energy con- 
strained by heavy-ion experiments; (f ) relationship between the transition density and the 
size of neutron-skin in 208 Pb; (g) systematics and constraints on the transition pressure at the 
core-crust boundary We will then study several global properties of neutron stars including 
the mass, radius, and the moment of inertia as well as their crustal fractions. The focus will 
be on effects of the density dependence of the symmetry energy on these observables. We 
will also check the inner crust EOS dependence of properties of neutron stars. 



5.1. Instabilities in neutron-rich matter and the core-crust transition density 

in neutron stars at /^-equilibrium 

In the subsection 12 . 31 we studied analytically the relationship between the dynamical and 
thermodynamical methods. To appreciate the relationship more clearly and quantitatively, 
we present here a numerical example in the p p vs. p n plane. First, it is important to recognize 
that the right-hand side of Eq. (13"2"|) just determines the thermodynamical instability of 
asymmetric nuclear matter. Shown in Fig. [3] are the boundaries of the mechanical (also 
known as the isothermal spinodal (ITS)), chemical (also known as the diffusive spinodal 
(DS)) and total instabilities without requiring the /3-equilibrium and charge neutrality using 
the MDI interac tion with x = a t zero temperature. A similar figure has been shown in our 
previous paper (jXu et al. 2008a ) but in the p ~ 5 or P ~ p plane at finite temperatures. 



Inside the ITS curve the system is mechanically unstable, while between the ITS curve and 
the DS curve the system is chemically unstable. The total instability is identified by the 
condition 



\ 2 

8p n dp p ( dp n 



dp n dp p \ dp 



< 0. (85) 
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Fig. 3. — (Color online) The mechanical, chemical and total instability boundaries shown in the 
pn ~ Pp plane using the MDI interaction with x = at zero temperature. The p n vs. p p for the 
npe matter at /^-equilibrium is shown as the dash-dot line. The core-crust transition density from 
the thermodynamical method is indicated with the filled dot. 

It is seen clearly that the total instability region obtained using the above condition covers 
the region of both mechanical and chemical insta bilities. This observation is co nsistent with 
the earlier finding by Margueron and Chomaz fjMargueron fc Chomazll2003h . In the npe 
matter when the /3-equilibrium and charge neutrality conditions are imposed, the p n and p p 
are correlated with each other. For the MDI interaction, this correlation can be obtained 
from the x p versus p curves shown in the window (c) of Fig. [2J With x — 0, this correlation 
is shown with the dash-dot line in Fig. [3j The cross point of this line and the boundary of 
total instability corresponds to the core-crust transition density in the npe matter within 
the thermodynamical approach. The density gradient term and the Coulomb interaction 
generally reduce slightly the instability region, thus the dynamical method normally leads 
to a slightly lower transition density. 



stable 
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It is necessary to note here that the onset of instab ilities has been associated with the 
so-called liquid-gas phase transition in nuclear matter ([Siemens I Il983l ). An experimental 
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manifestation of the liquid-gas phase transition is the well-known multifragmentation phe- 
nomenon in heavy-ion collisions, see, e.g., refs. fjChomaz et al.ll2004l ; iDas et al.ll2005l ) for a 
recent review. It is seen from Fig. [3] that both the dynamical and thermodynamical mod- 
els give a liquid-gas phase transition density of about 0.63p for symmetric (i.e., p n = p p ) 
nuclear matter (SNM) at T = 0. While the latter depends slightly on the interaction used, 
it is consistent with previous calculations on the boundary of mechanical instability in cold 
SNM, see, e.g., refs. fTMuller fc Serot lll995l : Ei fc Ko Ill997h . We thus conclude that both the 
dynamical and thermodynamical methods give the right asymptotical value for the transition 
density when one goes from the npe matter to the symmetric nuclear matter at T=0. 



5.2. Constraining the core-crust transition density in neutron stars 

We now turn to the numerical calculations and comparisons of the core-crust transition 
densities within both the dynamical and thermodynamical methods using the full EOS and 
its PA with the MDI and Skyrme interactions. The transition density can be directly ob- 
tained by carrying out the analyses in the p vs. x p plane. We stress again that in principle 
the transition density should be calculated from Eq. (fT2]) or Eq. (T2"T|) . and the /3-equilibrium 
condition should be expressed as Eq. ( l24|) . While in practical calculations, the parabolic ap- 
proximation has often been used in determining the /3-equilibrium condition using Eq. f[5"41) 
and in evaluating the V t her using Eq. ( 1231) . To first evaluate effects of the PA, we show in 
Fig. H] the density dependence of Vdyn and V[ her using the MDI interaction with x = and 
the Skyrme force R CT within both the dynamical and thermodynamical methods with the full 
EOS and its PA. Here, we have defined 

V ^r = V ther -^/[p—) (86) 

and it should be noted that V t ' her has the same vanishing point as the Vther and the same 
dimension as the Vd yn - For the MDI interaction with x = the transition densities using 
the full EOS within the dynamical and thermodynamical method are 0.065 fm -3 and 0.073 
fm -3 , respectively. While the corresponding results using the PA are 0.080 fm -3 and 0.090 
fm -3 , respectively. For the Skyrme force R CT the transition densities are 0.057 fm -3 and 0.066 
fm -3 using the full EOS, while the corresponding values with the PA are 0.084 fm -3 and 
0.093 fm~ 3 , by using the dynamical and thermodynamical method, respectively. Thus, the 
transition densities are generally lower with the dynamical method as we mentioned earlier, 
as the density gradient term and the Coulomb interaction make the system more stable. 
However, the PA significantly lifts the transition density regardless of the approach used. In 
fact, the difference between calculations using the full EOS and its PA is much larger than 
that caused by using the two different methods. 
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Fig. 4. — (Color online) The density dependence of V for the MDI interaction with x = (left 
window) and the Skyrme force Ro- (right window) using both the dynamical and thermodynamical 
methods with the full EOS and its parabolic approximation (PA). 



As we have mentioned in the introduction, it is well known that the tr ansition density de- 
pends sensitively on the E. qym (p) . Many interesting studies, see, e.g., refs. fjOyamatsu fc Iida 



2007I ; lLattimer fc Prakash ll2007l ; lDouchin fc Haensel II2000I ). have been carried out using var 



ious E sym (p). In the following, we present and compare the systematics of the transition 
density using the MDI interaction with the varying x parameter and 51 Skyrme forces. Since 
the density dependence of the symmetry energy can be well characterized by the L and K sym 
parameters, we examine the pt as a function of L and K sym in Fig. [5j Shown in the left 
panels are the pt as a function of L by using both the dynamical and thermodynamical 
methods with the full EOS and its PA. The same quantities are shown as a function of K sym 
in the right panels. It is interesting to see that both the dynamical and thermodynamical 
methods give very similar results with the former giving slightly smaller p t than the later 
(the difference is actually less than 0.01 fm -3 ) and this is due to the fact that the former 
includes the density gradient and Coulomb terms which make the system more stable and 
lower the transition density. The small difference between the two methods implies that the 
effects of density gradient terms and Coulomb term are unimportant in determining the pt- 
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Fig. 5. — (Color online) The transition density as a function of L (left panel) and K sym (right 
panel) by using both the dynamical and thermodynamical methods with the full EOS and its 
parabolic approximation. The MDI (upper windows) and Skyrme interactions (lower windows) are 
used. 



In addition, it is also interesting to see that the transition density decreases almost linearly 
with the increasing L especially in the calculations wit h the full EOS. This obs ervation is 
consistent with that found recently by Oyamatsu et al. (jOyamatsu Sz Iida 1 120071 ). We note 
here that there are some interactions with larger x values in the MDI interaction giving 
negative and/or very small values for the L parameter. These interactions with negative 
and/or very small values for the L parameter, however, lead to neutron-s kins in 208 Pb in- 
consistent with the existing data (jSteiner fc Li Il2005bl ; IChen et al. Il2005bl ). Since they are 
still somewhat theoretically interesting, we have thus also examined the possible transition 
density with these interactions. We find that for the interaction parameters with L < 7 MeV 
in the MDI interaction (x > 1.17), the transition density does not exist and the npe matter 
is always unstable. This is due to the fact that the symmetry energy is so soft that the ^ L 



is always negative while the ^ is always positive at low densities, and thus the stability 

app 



condition f^§±£ 

Opn ap p 



dp n 

dp p 



> can never be satisfied. 
It is clear from all existing calculations that the p t is sensitive to the density dependence 
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of the nuclear symmetry energy. The latter can be well characterized by the slope L and the 
curvature K sym . Naturally, there are some correlations between the L and K sym determined 
by the interaction energy density functional used. For the MDI interaction, the L and K sym 
both change linearly with the parameter x. Therefore they are linearly correlated. Similarly, 
the L and K sym also correlated within the SHF model. It is therefore not surprising that 
the variation of pt with K sym is very similar to that with L, as shown in the right panels of 
Fig. El 

We now apply the experimentally constrained L to the p t — L correlation obtained 
using the full EOS within the dynamical method in constraining the pt- In our earlier 



transport model studies of the isospin diffusion data in heavy-ion reactions (|Chen_et_al 



2005al ; iLi &: Chen 1 120051 ) . the complete MDI interaction was used. The extracted L value is 



88 ± 25 MeV if one defines the E sym using the PA in Eq. ([35]) or 86 ± 25 MeV if one uses 
the exact expression of ^^(Eq. (1491) ) corresponding to the full MDI EOS. Using the latter 
in comparison with the full MDI results shown in Fig. [6l we conclude that the transition 
density is between 0.040 fm -3 and 0.065 fm~ 3 . This constrained range is significantly below 
the fiducial value of p t ~ 0.08fm~ 3 often used in the literature and the estimate of 0.5 < 



Pt/ Po < 0.7 made in ref. (ILattimer &; Prakash 1120071 ) within the thermodynamical approach 
using the parabolic approximation of the EOS. This difference is understandable as we shall 
explain in detail below. 



5.3. Understanding effects of the parabolic approximation of the EOS on the 

core-crust transition density in neutron stars 

It is also seen from Fig. \5\ that except at very small values of L and K sym , there is a 
big difference between results obtained using the full EOS and its parabolic approximation 
within both the dynamical and thermodynamical methods. Especially at high values of 
L and K sym , the p t from the PA increases while t he one from t he full EOS continuously 



decreases. We also notice that Kubis's calculations (IKubis 1 120071 ) coincide with our results 
using the MDI interaction within the thermodynamical method with the PA. Why is the 
transition density so sensitive to whether one used the PA or not? To answer this question, 
we first notice that both the first and second derivatives of the EOS are involved in the 
stability conditions. The EOS can be expanded according to E b (p, x p ) up to the fourth order 
term of (1 — 2x p ) according to 

E b (p,x p ) = E (p) + E sym (p)(l - 2x p f + E 4 (p)(l - 2x p f 

+ 0(1 -2x p f. (87) 
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Fig. 6. — (Color online) The relation between the transition density and L by using the thermody- 
namical method with the MDI (left panel) and Skyrme (right panel) interactions. For both kinds 
of interactions, up to the 2nd and 4th orders in isospin asymmetry are used in the expansions of 
the corresponding full EOS. 



Only even order terms of (1— 2x p ) appear as the strong interaction is assumed to be symmetric 
for exchanging neutrons with protons. The first and second order derivatives of the energy 
with respect to respectively, 

OE 

= -AE sym {p){l - 2x p ) - 8£ 4 (p)(l - 2x p f + 0(1 - 2x p f . (88) 

^ = 8E syrn (p) + 48£ 4 (p)(l - 2x p f + 0(1 - 2x p f. (89) 

At /3-equilibrium the npe matter is usually highly neutron rich, so the (1 — 2x p ) is not far 
from 1. Thus the higher order terms in (1 — 2x p ) are normally not negligible. Moreover, 
although the coefficient E4 is normally smaller than the E sym , the contribution to the ^ k 

and the from the E4 term gains a factor of 2 and 6, respectively, compared to that from 
the E sym term. Thus, mathematically one expects the £4 term to play an important role in 
locating the stability boundaries in asymmetric nuclear matter and the core-crust transition 
density in neutron stars. It is also easy to understand why the effect is stronger with stiffer 
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symmetry energy functionals. At sub-saturation densities near the pt, the proton fraction 
x p is lower with the stiffer symmetry energy. It is the opposite at supra-saturation densities. 
A numerical example can be found in the window (c) of Fig. [2] for the MDI interaction with 
x = (softer) and x — — 1 (stiffer). It is seen that with the stiffer symmetry of x = —1, 
near the transition density (1 — 2x p ) is indeed larger than that with a softer symmetry 
energy of x = 0. Therefore, a larger error will be introduced in calculating the pt using the 
parabolic approximation with stiffer symmetry energy functionals. To be more quantitative, 
we compare in Fig. O the p t as a function of L obtained within the thermodynamical method 
using the full EOS with those obtained by expanding the EOS to the second and 4th orders 
in (1 — 2x p ). The left window is for calculations with the MDI interaction, and the right 
one with the 51 Skyrme forces. We notice that the convergence speed is very slow, and 
not only the fourth order term but also the sixth, eighth or even further higher order terms 
should be considered (For the Skyrme forces, we noted the calculations up to the 8th order 
approximation still leads to a significant error compared to the full EOS). We also notice here 
that the EOS of asymmetric nuclear matter always contains the higher-order terms in isospin 
asymmetry at least due to the kinetic contribution. Moreover, if we use the Eq. (135]) instead 
of Eq. (I3T)|) in calculating the symmetry energy and then reconstruct the EOS as E(p, 5) = 
E(p, 5 = 0) + [E(p, 5=1) —E(p, 5 = 0)]5 2 + 0(5 A ) as in the parabolic approximation, almost 
the same transition densities are obtained as the second order approximation shown in Fig. 
[6j Our results thus indicate that one may introduce a huge error by assuming a priori that 
the EOS is parabolic for a given interaction in calculating the p t . It is thus clear that the 
correct transition density can hardly be obtained without knowing the exact expression of 
Ej,(p, x p ) for a given interaction. Interestingly, these features agree with the early finding 
(lArpoen Ill972l ) that the p t is very sensitive to the fine details of the nuclear EOS. 



5.4. Correlation between the core-crust transition density in neutron stars 

and the size of neutron-skin in 208 Pb 



It is also well known that the sizes of neutron skins in heavy nuclei are sensitive to the 
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). However, available data of neutron-skin thickness obtained us- 



ing hadronic probes are not accurate enough yet to constrain tightly the symmetry energy. 
Interestingly, the parity radius experiment (PREX) at the Jefferson Laboratory aiming to 
measure the neutron radius o f 208 Pb via parity viol ating electron scattering (Jefferson Labo- 
ratory Experiment E-00-003) ([Horowitz et al.ll200ll ) hopefully will provide much more precise 
data in the near future. It can then potentially constrain the symmetry energy at low densi- 
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Fig. 7. — (Color online) The slope parameter L and the transition density pt as a function of the 
neutron skin thickness S of 208 Pb using the dynamical method with the full EOS with the Skyrme 
interactions. The solid lines indicate the linear fits. 



ties and thus the core-crust transition density more stringe ntly. It has been shown by many 
authors and also in our previous work (IChen et al. Il2005bl ) that the neutron skin thickness 
S increases linearly with L. Given the fact that the transition density p t decreases almost 
linearly with the increasing L as shown above, it is interesting to e xamine the correlation be- 



tween the S and p t . Such kind of study was first carried out in ref. (IHorowitz &: Piekarewicz 



20011 ) using the RMF EOS and the p t calculated within the RPA approach. 



Shown in Fig. [7] are the p t and L versus the neutron skin thickness S of 208 Pb obtained 
by using the dynamical method with the full SHF E OS. As known befo re, the neutron 
skin thickness increases linearly with the increasing L ( IChen et al. Il2005bl ). Moreover, the 
transition density shows a decreasing trend with the increasing neutron skin thickness. As 
the p and E sym (p ) are different for the various sets of Skyrme forces, the data points do not 
show a very strong linear correlation. However, the tendency is clea r. This trend is consistent 



with the RPA calculations u sing the RMF EOS's by Horowitz et al. (IHorowitz fc Piekarewicz 



200ll ; ICarriere et al. 1120031 ). 
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Fig. 8. — (Color online) The transition pressure Pt as a function of pt and L within the ther- 
modynamical method with the full EOS and its parabolic approximation using the MDI (upper 
windows) and Skyrme (lower windows) interactions. 



5.5. The pressure at the inner edge of neutron star crust 



The pressure at the inner edge is an important quantity related directly with the crustal 
fraction of the moment of inertia which ca n be measurable indirectly from observations of 
pulsar glitches (ILattimer &: Prakash 1120071 ). In principle, having determined the transition 
density it is straightforward to calculate the corresponding pressure using the formalisms 
outlined in the subsection 13.31 Before presenting the numerical results , it is very instructive 



to qu ote the analytical estimation obtained by Lattimer and Prakash (ILattimer fc Prakash 
20071 ) for the transition pressure 



Pi 



I<Qf>l (fk_ 



9 Po VPo 



pA 



I -St 



E S ym(pt) 



dE sym {p) 

dp 



pt 



(90) 



where Kq is the incompressibility of SNM at po and St is the isospin asymmetry at pt- Besides 
the implicit dependence on the symmetry energy through the pt and St, the Pt also depends 
explicitly on the value and slope of the E sym (p) at p t . Thus the P t depends very sensitively 
on the E sym (p). Noticing that the Eq. (I90p is an estimate using the thermodynamical method 
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Fig. 9. — (Color online) Same as Fig. [8] but within the dynamical method. 

with the PA, it is interesting to compare numerically its predictions with the results obtained 
according to the original formalisms in the subsection I3.3I using both the full EOS and its 
PA within the thermodynamical and dynamical method, separately. 

In Fig. [8] we show the P t as a function of pt (left windows) and L (right windows) by 
using the thermodynamical method with and without the parabolic approximation. The 
same quantities with the dynamical method are shown in Fig. [9j Both the MDI (upper 
windows) and Skyrme (lower windows) interactions are used. The results from Eq. (190]) using 
the p t and E sym corresponding to the full EOS and its PA are also shown for comparisons. 
It is interesting to see that the Eq. (19"U|) predicts qualitatively the same but quantitatively 
slightly higher values compared to the original expressions for the pressure with or without 
the PA for both the thermodynamical and dynamical methods even though this formula was 
derived from the thermodynamical method using the PA. This observation is consistent with 
the results shown in the window (b) of Fig. [2l namely, the direct effect of using the full EOS 
or its PA on the pressure is small although the PA may affect strongly the transition pressure 
P t by changing the transition density p t . The P t essentially increases with the increasing p t 
in calculations using the full EOS, but a complex relation between the P t and p t is obtained 
using the PA. The observed large difference in P t is due to the strong PA effect on the p t . 
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Fig. 10. — (Color online) Pt as a function of pt by using dynamical method without parabolic 
approximation for both MDI interaction and SHF calculations. The shaded band represent 
the constraint from the isospin diffusion data. 



Moreover, the latter does not vary monotonically with L for the PA as shown in Fig. [5j Thus 
the PA of the EOS leads to a very different P t compared to the calculations with the full 
EOS especially for the stiffer symmetry energy functionals. 

It is also interesting to examine the range of Pt corresponding to the pt and L constrained 
by the heavy-ion reaction data. In Fig. [10l we show the Pt as a function of pt by using the 
dynamical method and the full EOS for both the MDI (solid line) and the Skyrme (filled 
squares) calculations. It is interesting to see that the MDI and Skyrme interactions give 
generally quite consistent results. Corresponding to the pt constrained in between 0.040 
fm~ 3 and 0.065 fm -3 , the P t is limited between 0.01 MeV/fm 3 and 0.26 MeV/fm 3 with the 
MDI interaction as indicated by the shaded area, which is s ignificantly less than the fid ucial 



value of Pt ~ 0.65 MeV/fm often used in t he literature (ILattimer fc Prakash 1 120071 ). As 



pointed out in a recent work by Avancini et al f Avancini et al. 2008bl ). the value of Pt ~ 0.65 



MeV/fm 3 may be too large for most mean-field calculations without the PA. We notice here 
that among the 51 Skyrme interactions listed in Tables [1] and [2J the following 7 interactions, 
i.e., the SkMP, SKO, R CT , G a , SkI2, SkI3, and SkI5, are consistent with the constraints from 
heavy-ion reactions. 
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Fig. 11. — (Color online) pt and Pt as a function of L by using dynamical method for MDI 
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In closing this subsection, we examine how the transition density and pressure may be 
sensitive to variations of the coefficients D pp = D nn = D np in the MDI interaction. As we 
have pointed out in subsection l2.ll for the dynamical method, we introduced phenomenolog- 
ically the empirical values of D pp = D nn = D np = 132 MeV-fm 5 for the gradient coefficients 
in the MDI interaction. These values are obviously not obtained self-consistently. Shown in 
Fig. [TT] are the pt and Pt as functions of L by using the dynamical method with the full MDI 
EOS but different values of the coefficients D pp = D nn = D np , namely, D pp = D rm = D np = 
50, 132, and 200 MeV-fm 5 , respectively. We note from Fig. [TT] that changing the value of 
D pp = D nn = D np from 50 to 200 MeV-fm 5 leads to at most a variation of about 0.007 fm~ 3 
for p t and 0.06 MeV/fm 3 for P t . These results thus indicate that the transition density and 
pressure are rather insensitive to the variation of D pp = D nn and D np = D pn . 



5.6. Constructing the EOS from the center to the surface of neutron stars 

With a clear understanding about the core-crust transition density as we discussed 
above, we now investigate several other properties of the crust and the whole neutron star. 
To proceed, it is necessary to know the EOS over a broad density range from the center 
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to the surface of neutron stars. Besides the possible appearance of nuclear pasta in the 
crust, various phase transitions and non-nucleonic degrees of freedom may appear in the 
core. In this work, we restrict ourselves to the simplest and traditional model. We make 
the minimum assumption that the core contains the uniform npefi matter only and there 
is no phase transition. Results of this kind of calculations serve as a useful baseline for 
understanding general features of astrophysical observations. Significant deviations from 
observations may indicate the onset of non-traditional physics. 

For the core we use the MDI EOS and its PA shown in Fig. [3 In the inner crust 
of densities between p out and p t where the nuclear pastas may exist, because of our poor 
knowledge about its EOS from first principle, following Carriere et al. (jCarriere et al. 1120031 ) 
we construct its EOS according to 

P = a + b e 4/3 . (91) 

This polytropic form with an index of 4/3 has been found to be a good approximation to 
the crust EOS fjLink et al. Ill999l ; lLattimer fc Prakash 1 120001 ) and we will discuss how our 
results are sensitive to the polytropic index later. The p out = 2.46 x 10~ 4 fm -3 is the density 
separating the inner from the outer crust. The constant a and b are determined by 

p , 4 / 3 _ p J/ 3 
r out t t r t t out 



4/3 4/3 



-out 



Pt-P, 



Oil t 



4/3 4/3 



(92) 



— e 



out 



where P t , e t and P ou t, e out are the pressure and energy density at p t and p„ n t, respectively. 
In the outer crust with 6. 93 x 1CT 13 fm~ 3 < p < p out we use the EOS of BPS (IBaym et al. 
1971al : llida &: Satol ll997h , and in the region of 4.73 x 10" 15 fm" 3 < p <6.93 x 10" 13 fm" 3 



we use the EOS of FMT f lBavm et al. Ill971al ). 



Shown in Fig. [12] are the selected EOS for the different parts of the neutron star. 
As we have discussed earlier, the pt is obtained by studying the onset of instabilities in 
the core, namely it is the critical density below which small density fluctuations will grow 
exponentially. The p t is thus determined by the EOS of the core only. We use here the p t 
obtained within the dynamical method using the full EOS and its parabolic approximation 
with the MDI interaction of x = and x = — 1. The corresponding values of et are indicated 
by the vertical lines in Fig. [12j Using the above combination of EOS's for the different parts 
of the neutron star, the radial distribution of the total energy density and the pressure in 
neutron stars is continuous as required, but the derivative of the pressure is not continuous 
at p t and p ou t- It is seen that the EOS for the inner crust is quite different using the Full 
EOS or its PA especially with x = —1. Interestingly, one can see that the famous BPS EOS 
extended to the inner crust is between the parameterized EOS's with x = and x — — 1. 
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Fig. 12. — (Color online) The EOS of different parts of neutron stars. The energy density at pt 
and pout is indicated in the figure as et and e ou t, respectively, and the full and PA results of MDI 
interaction with x = and x = — 1 are shown. 



5.7. The mass-radius correlation of neutron stars 



With the EOS constructed above, in the next three subsections we study several key 
properties of the crust and the whole neutron star using the formalisms outline in section HI 
We carry out numerical calculations for all interested quantities. For the crustal fraction 
of the moment of inertia, we also compare our numerical calc ulations with predictions of 
the a nalytical expression put forward by Lattimer and Prakash (ILattimer &: Prakash 1 120071 . 
2000|) . In this subsection, we focus on effects of the symmetry energy on the mass-radius 
correlation. We use the MDI interaction with x = and x = —1 consistent with the existing 



heavy- ion reaction data (ILi et al. 1120081 ) 



The resulting mass-radius correlation is shown in Fig. [13j For the softer symmetry 
energy (x = 0) the M decreases with increasing R, while for the stiffer symmetry energy(x = 
— 1) the radius remains almost unchanged or even decreases with decreasing mass near 
R = 13.5 km. For M > 0.53M Q the radius is larger for x = —1, while for M < O.53M 
the radius is larger for x = 0. For nucleonic matter, a stiffer symmetry energy leads to a 
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Fig. 13. — (Color online) The M-R relation of static neutron stars from the full EOS and 
its parabolic approximation as well as the test case (see text for details) with the MDI 
interaction with x = and x = —1. For the Vela pulsar, the constraint of AI/I > 0.014 
implies that allowed masses and radii lie to the right of the line linked with solid squares 
(p t = 0.065 fm -3 and P t = 0.26 MeV/fm 3 , obtained i n the present wor k) or open squares 



(p t = 0.075 fm" 3 and P t = 0.65 MeV/fm 3 , used in ref. flLink et al. Ill999h ). 



stiffer EOS for the liquid core, but a lower core-crust transition density. The crossing point 
of the M-R curves with x = and x = — 1 is a result of this competition. It is clearly shown 
that with the PA the radius is larger at a fixed mass especially for the stiffer symmetry 
energy of x = — 1. To better understand the role of the transition density in determining 
the M-R relation, we also made an additional test by using the full MDI EOS but with the 
p t obtained from using the PA. The results are shown with the dotted lines. They are very 
close to the results obtained consistently using the PA in calculating both the EOS and the 
transition density. Thus, the mass-radius relation, especially the radius, seems to be quite 
sensitive to the location of the inner edge. The small difference between the full EOS and 
its PA for the core (shown in the panels (a) and (b) of Fig. [2J) has a negligible effect on the 
M-R relationship once the inner edge is fixed. These features can be seen more clearly in 
Fig. E] where the mass and radius are displayed separately as functions of the central energy 
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Fig. 14. — (Color online) The mass and radius of neutron stars as functions of the central energy 
density using the MDI interaction with x = and x = — 1. The results from three methods are 
shown for comparison. See text for details. 



density. Very similar masses are obtained independently of how the pt was calculated for 
a given x parameter. However, the radii are appreciably different for the stiffer symmetry 
energy with x = — 1 using the full EOS or its PA because of their different p t values. Also, 
since the test case has the same p t as the PA, it thus leads to the same radii as the PA for 
both x = and x — — 1. 



5.8. The crust thickness and the core size of neutron stars 



For a given neutron star of total mass M and radius R, what are the respective sizes 
of its core and crust? How do they depend on the stiffness of the symmetry energy? How 
do they depend on the neutron star mass Ml It is well known that the size of neutron skin 
in heavy nuclei increases with the increasing L as shown in Fig. [7J How does the thickness 
of neutron star crusts depend on the LI These are among the interesting questions we shall 
discuss in this subsection. First, we display in Fig. [15] the radial energy density profile for 
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Fig. 15. — (Color online) The radial energy density distribution of neutron stars, using MDI 
interaction with x = and x = — 1, at total mass 0AM@, 1.0-Mq and IAMq. The interface 
between uniform part and crust part is indicated. R t is the radius of the liquid core and et is the 
energy density at the edge of the liquid core and the crust. 



neutron stars of total mass O.4M , 1.0 Mq and 1.4M using the MDI interaction with x = 
and x — — 1, respectively. The inner edge separating the uniform core from the crust is 
indicated by the vertically dotted lines. The corresponding energy density e t is shown as the 
longitudinally dotted lines. It is very interesting to see that the radius of the core increases 
while the thickness of the crust decreases with the increasing neutron star mass M. The 
lighter neutron stars generally have thicker and more diffusive crusts due to the competition 
between the gravitation and the nuclear forces. Moreover, this feature is independent of the 
symmetry energy used. It is also seen that the stiffer symmetry energy with x = —1 predicts 
a larger core but a thinner crust for a given mass M. More quantitatively, for a canonical 
neutron star of M = 1.4M , the radius of the core is 10.89 km with x = and 12.55 km 
with x = —1, and the thickness of the crust is 1.09 km with x = and 0.72 km with x = — 1, 
respectively. Therefore, with a softer symmetry energy, a light neutron star can have a big 
radius due to its very thick crust. 
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Fig. 16. — (Color online) The whole radius R, the crust thickness AR, the core radius Rt as 
functions of L at fixed total mass of 0.4M©, I.OMq and IAMq, respectively. 



To study more systematically effects of the symmetry energy, we show in Fig. [16] the 
core radius Rt, the crust thickness AR and the total radius R as functions of L for a fixed 
total mass of O.4M , I.OMq and IAMq, respectively. It is seen that the R t increases almost 
linearly with the increasing L. The R t also increases with the increasing mass at a fixed L. 
This is because the stiffer the symmetry energy is, the larger the contribution of the isospin 
asymmetric part of the pressure will be, which makes the Rt larger. Moreover, the AR 
decreases with the increasing L especially for light neutron stars, as the transition density 
decreases with the increasing L. As the thickness of the crust AR and the core radius 
Rt depend oppositely on L, the total radius R = R t + AR of the neutron star may show 
a complicated dependence on L. We stress here that this is the result of a competition 
between the repulsive nuclear pressure dominated by the symmetry energy contribution and 
the gravitational binding. Interestingly, it is often mentioned that the crust of neutron stars 
bears some analogy with the neutron-skin of heavy nuclei. However, they show completely 
opposite dependences on the L. Namely, the size of neutron-skin usually increases with 
the increasing L as a result of the competition between the nuclear surface tension and 
the pressure difference of neutrons and protons, while the thickness of neutron star crusts 
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Fig. 17. — (Color online) The crustal fraction of neutron mass AM/M, the moment of inertia I of 
the whole star and the crust contribution AI as a function of L, at fixed total mass OAMq, I.OMq 
and IAMq, respectively. 



decreases with the increasing L as a result of the competition between the nuclear pressure 
and the gravitational binding. 



5.9. The crustal fractions of neutron star masses and moments of inertia 

What is the crustal fraction AM/M of the total mass and how does it depend on the 
symmetry energy? Since the mass is simply the integration of the energy density, from 
the profile of the energy density shown in Fig. H51 we expect the AM/M and AR/R have 
very similar dependences on L. Shown in the right window of Fig. [17] is the AM/M. The 
fractional mass of the crust decreases with the increasing L at a fixed total mass, and it 
decreases with the increasing total mass at a fixed value of L. The moment of inertia is 
determined by the distribution of the energy density. From the middle window, it is seen 
that the total moment of inertia increases with the increasing mass at a fixed value of L 
and increases with the increasing L at a fixed total mass. The dependence on L is relatively 
weak especially for the light neutron stars. However, the crust contribution of the moment 
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Fig. 18. — (Color online) The relation between the crustal fraction of the moment of inertia and 
the total mass or the radius of neutron stars, using MDI interaction with x = and x = — 1. The 
results of PA or exact method from direct calculation or LP's formula are shown for comparison, 
and the constraint of AI/I is also indicated. 



of inertia varies much more quickly with L. It decreases with the increasing neutron star 
mass at a fixed value of L and decreases with the increasing L at a fixed total mass. These 
are all consistent with the behaviors of the fractional mass and size of the crust. 

The crustal fraction of the moment of inertia AI/I is particularly interesting as it can 
be inferred from observations of pulsar glitches, the occasional disruptions of the otherwise 
extremely regular pulsations from magnetized, rota ting neutron stars. It can be expressed 
approximately as (ILattimer fc Prakash 1 120071 . l2000l ) 



AJ 



2SnP t R 3 (1 - 1.67f - 0.6£ 2 



3Mc 2 



i + 



2P f (l + 5£-14£ 2 
p t mc 2 ^ 2 



(93) 



where m is the mass of baryons and £ = GM/Rc 2 . A numerical verification of this formula 
is useful. Predictions of this formula (thin lines) are compared in Fig. [TS] with our direct 
numerical calculations (thick lines). Very interestingly, the analytical formula reproduces 
very well our results from direct numerical calculations using both the full EOS and its PA. 
Comparing calculations using the full EOS and its PA, one sees clearly big differences, again 
due to the corresponding differences in the transition density. For instance, using either the 
direct numerical calculation or the formula (1931) . at a fixed total mass M the AI/I increases 
using the full EOS while it decreases using the PA when the x parameter is changed from 
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x = — 1 to x = 0. As it was stressed in ref. (ILattimer fc Prakash Il2000l ). the AI/I depends 
sensitively on the symmetry energy at sub-saturation densities through the P t and p t , but 
there is no explicit dependence upon the EOS at higher-densities. 

Experimentally, the crustal fraction of the moment of inertia has been constrained as 
AI/I > 0.014 from studying the glitches of the Vela pulsar (ILink et al. 1119991 1 . As indicated 
in Fig. [TBI this limits the masses and radii of the neutron star. For example, from Fig. [18j it 
is indicated that the maximum mass is about 1.57M (0.73M Q ) while its minimum radius is 
about 11.6 (13.4) km for the MDI interaction with x = (x = —1) if the dynamical method 
is used to determine the Pt and pt- We note here that the very small mass for Vela pulsar 
constrained by this condition using the MDI with x — —1 is due to the associated small 
transition density and pressure. Combining the observational constraint of AI/I > 0.014 
with the upper bounds of p t = 0.065 fm -3 and P t = 0.26 MeV/fm 3 inferred from heavy-ion 
reactions, we can obtain a minimum radius of R > 4.7 + 4.0M/M Q km for the Vela pulsar. 
This limit is indicated by the solid squares in Fig. [131 According to this constraint, the 
radius of the Vela pulsar is predicted to exceed 10.5 km should it have a mass of 1.4M . It is 
worth mentioning that a constraint of R > 3.6 + 3.9.U/.U. km for this pulsar (see the open 
squares in Fig. [T3"j) has been derived previously in ref. (ILink et al. Ill999l ) by using p t = 0.075 
fm -3 and P t = 0.65 MeV/fm 3 . The difference between this and our prediction is due to the 
different p t and P t . 



5.10. The inner crust EOS dependence of neutron star properties 



As discussed in Eq. (|9"T|) of subsection 15.61 we have adopted the polytropic EOS of P = 
a + be 1 with 7 = 4/3 for the inner crust in the above calculations. This particula r polytropic 



EOS has been extensively used for studying the inne r crust in the literature (ILink et al 



19991 ; ILattimer fc Prakash 1 12000| ; ICarriere et al. 1 120031). However, due to the complexity of 



the inner crust, its EOS is rather uncertain (INegele fc Vautherin Ill973l ). Thus, it would be 
interesting to investigate how our results may be sensitive to the polytropic index 7. 

Firstly, let us see how the polytropic index affects the mass and radius of a neutron 
star. Shown in Fig. [19] is the M-R relation obtained using the full MDI EOS with x = 
and x = —1 with 3 different values of the polytropic index 7, i.e., 7 = 1.5,4/3, and 1. It is 
seen that the polytropic index has very little effects on the mass of neutron stars. On the 
other hand, it is interesting to see that the neutron star radius increases significantly with 
the deceasing polytropic index 7 especially for the softer symmetry energy (x — 0). This is 
due to the change of the crust thickness from varying the inner crust EOS. The observed 
symmetry energy dependence of the polytropic index effects on the neutron star radius can 
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Fig. 19. — (Color online) The M-R relation of static neutron stars from the full MDI EOS 
with x = and x = — 1. The different values of the polytropic index 7 in P = a + 6e 7 for 
the inner crust EOS, i.e., 7 = 1.5,4/3, and 1 have been used. 



be easily understood since the stiffer symmetry energy leads to a thinner thickness of the 
crust as shown in Fig. [16] and thus less sensitivity to the variation of the inner crust EOS. 
Our results thus indicate that for softer symmetry energies, an accurate inner crust EOS is 
important for the precise determination of the neutron star radius. 

In order to see the inner crust EOS dependence of the crustal fraction of the moment 
of inertia A///, we show in Fig. [20] the AI/I as functions of the total mass and the radius 
of neutron stars from the full MDI EOS with x = and x = — 1 using 7 = 1.5,4/3, and 1. 
Noting the very weak dependence of the neutron star mass on the 7 index, we can see clearly 
from the left window of Fig. [20] that the AI/I is not so sensitive to the variation of the inner 
crust EOS, especially for stiffer symmetry energies. We notice here that the pronounced 7 
index dependence of AI/I as a function of the neutron star radius is due to the fact that 
the neutron star radius depends significantly on the 7 index as shown in Fig. [TUl especially 
for the softer symmetry energy {x = 0). 



Al so indicated in Fig . [20l is the constraint of AI/I from studying the glitches of the Vela 
pulsar (ILink et al. Ill999l ). It is very interesting to see that the neutron star mass and radius 
at AI/I = 0.014 exhibit a very weak dependence on the polytropic index 7. This nice feature 
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Fig. 20. — (Color online) The relation between the crustal fraction of the moment of inertia and 
the total mass or the radius of neutron stars, using MDI interaction with x = and x = — 1. The 
different values of the polytropic index 7 in P = a + be 1 for the inner crust EOS, i.e., 7 = 1.5, 4/3, 
and 1 have been used. The constraint of AI/I is also indicated. 



implies that the obtained constraint on the minimum radius of R > 4.7 + 4.0M/M Q km for 
the Vela pulsar in the present work is not sensitive to the inner crust EOS and still holds. 
Moreover, we can see from Fig. 1201 that the robustness of the constraint R > 4.7 + 4.OM/M 
km against the variation of the inner crust EOS is actually due to the very small value of 
AI/I for the Vela pulsar. For higher values of AI/I, on the contrary, the constraint will 
depend significantly on the inner crust EOS. 



6. Summary 

In summary, we first analyzed the relationship between the well established dynamical 
and thermodynamical methods for locating the inner edge separating the uniform liquid core 
from the solid crust in neutron stars. It is shown analytically that the thermodynamical 
method corresponds to the long- wavelength limit of the dynamical one when the Coulomb 
interaction is neglected. Moreover, it is shown that the results obtained from using the full 
expression of the EOS for a given interaction are very similar for the two methods. However, 
the widely used parabolic approximation to the EOS of asymmetric nuclear matter leads 
systematically to significantly higher core-crust transition densities and pressures, especially 
for stiffer symmetry energy functionals regardless of the specific method used in calculating 
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the transition density. Our results thus indicate that one can hardly obtain the accurate 
transition density without knowing the complete EOS and may introduce a huge error by 
assuming a priori that the EOS is parabolic in isospin asymmetry for a given interaction. 
Based on systematical calculations using the modified Gogny force (MDI interaction) and 
selected 51 Skyrme interactions widely used in the literature, it is shown that the transition 
density and pressure are very sensitive to the density dependence of the nuclear symmetry 
energy. We also systematically investigated several properties of neutron star crust. We 
found that the thickness, fractional mass and moment of inertia of neutron star crust are 
all very sensitive to the slope parameter L of the nuclear symmetry energy through the 
transition density p t and the results depend on whether one uses the full EOS or its parabolic 
approximation. Therefore, accurate knowledge on the nuclear symmetry energy at sub- 
saturation densities is required to fully understand the properties of neutron star crusts. 

Using the MDI EOS of neutron-rich nuclear matter constrained by the recent isospin 
diffusion data from heavy-ion reactions in the same sub-saturation density range as the 
neutron star crust, the transition density and pressure at the inner edge of neutron star 
crusts are limited to 0.040 fm~ 3 < p t < 0.065 mi" 3 and 0.01 MeV/fm 3 < P t < 0.26 MeV/fm 3 , 
respectively. The constrained range of the transition density is significantly below the fiducial 
value of p t ~ 0.08 fm~ 3 often used in the literature and the estimate of 0.5 < pt/po < 0.7 
made previously within the thermodynamical approach using the parabolic approximation 
of the EOS while that of the Pt is also significantly less than the fiducial value of Pt ~ 0.65 
MeV/fm 3 often used in the literature. The newly constrained transition density and pressure 
together with the condition AI/I > 0.014 for the crustal fraction of the moment of inertia 
extracted from studying glitches of the Vela pulsar allow us to set a new limit on the radius 
of the Vela pulsar, i.e., R > 4.7 + 4.OM/M . It is significantly different from the previous 
estimate and thus puts a new constraint for the mass-radius relation of neutron stars. 

Finally, it is worth noting that in the present work, we have only considered the non- 
accreting crusts of cold, non-rotating nucleonic neutron stars. In the next step, we plan to 
extend the study to accreting neutron stars. It will be especially interesting to investigate 
how the finite temperature, the strong magnetic field and neutrino trapping may affect the 
transition density and pressure reported here. Moreover, there are still many interesting 
issues regarding the neutron star crust, such as its composition, thermal, transport and 
mechanical properties that are important for a better understanding of the structure and 
evolution of protoneutron stars, the x-ray bursts and the emission of gravitational waves from 
neutron stars. More information from terrestrial nuclear reactions especially those induced 
by radioactive beams will certainly contribute to resolving these issues. 
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